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We develop two numerical methods to solve the differential equations with deviating 
arguments for the motion of two charges in the action-at-a-distance electrodynamics. Our 
first method uses Sturmer's extrapolation formula and assumes that a step of integration 
can be taken as a step of light ladder, which limits its use to shallow energies. The 
second method is an improvement of pre-existing iterative schemes, designed for stronger 
convergence and can be used at high-energies. 
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1. Action-at-a-distance electrodynamics 

In this paper we continue an investigation initiated in and here extended to the 
numerical analysis of the 3-dimensional two-body problem in the relativistic action- 
at-a-distance electrodynamics of Wheeler and Feynmana , henceforth called 3D WF. 
The dynamics follows from the Schwarzschild-Tetrode-FokkerB direct-interaction 
functional. Equations of motion are derived from Hamilton's principle for the action 
integral 

S = - J2I m % cdsi - (eiZj/c) J J & [\\xi - Xj\\ ) X% ■ x ds t dsj, 

where the four- vector Xi(.Si) represents the four-position of particle i parametrized 
by arc- length s^, double bars indicate four- vector modulus \\xi — Xj\\ = {xi — Xj) ■ 
(xi — Xj) and the dot indicates Minkowski relativistic scalar product of four- vectors 
with metric tensor = diag{\, —1,-1,-1). Integration is to be carried over the 
whole particle trajectories, at least formally 0, unlike other usual additive cases 
of classical mechanics. The above action integral describes an interaction at the 
advanced and retarded light-cones with an electromagnetic potential giyen by half 
the sum of the advanced and retarded Lienard-Wierchert potentials Wheeler 
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and Feynman showed that electromagnetic phenomena can be described by this 
direct action-at-a-distance theory in complete agreement with Maxwell's theory as 
far as the classical experimental consequencesB'tl This direct-interaction formu- 
lation of electrodynamics was developed to avoid the complications of divergent 
self-interaction, as there is no self-interaction in this theory, and also to eliminate 
the infinite number of field degrees of freedom of Maxwell's theoryu. It was a jpeat 
inspiration of Wheeler and Feynman in 1945, that followed a lead of TetrodeH and 
showed that with the extra hypothesis that the electron interacts with a completely 
absorbing universe, the advanced response of this universe to the electron's retarded 
field arrives at the present time of the electron and is equivalent to the local instan- 
taneous self-interaction of the Lorentz-Dirac theoryu plus the retarded interaction 
among chargesEI The action-at-a-distance theory is also symmetric under time re- 
versal, as the Fokker action includes both advanced and retarded interactions. Dis- 
sipation in this time-reversible theory is due to interaction with the other charges of 
the universe and becomes a matter of statistical mechanics of absorptions. The area 
of Wheeler- Feynman electrodynamics has been progressing slowly but steadily since 
1945: quantization was achieved by use of the Feynman path integral technique and 
the effect of spontaneous emission was successfully described in terms of interaction 
with the future absorber, in agreement with quantum electrodynamics. It was also 
shown that it is possible to avoid the usual divergencies associatecL with quantum 
electrodynamics by use of proper cosmological boundary conditions!!! As far as un- 
derstanding of the dynamics governed by the equations of motion, the state of the 
art is as follows: the exact circular orbit solution to the attractive two-body prob- 
lem was proposed in 1946 by Schonberglill and rediscovered in 1962 by Schildta. 
The 1-dimensional symmetric two-electron scattering is a special case where the 
equations of motion simplify a Jotand it has been studied by many authors, both 
analytically and numerically E3oll3. In this very special case the initial value func- 
tional problem surprisingly requires much less than an arbitrary initial function to 
determine a solution manifold with the extra condition of bounded manifold for all 
times. It was shown that the solution is uniquely determined by the inter-electronic 
distance at the turning point if this distance is large enougnlil. As a result of this 
theorem, there is a single continuous parameter (the positive energy) describing the 
unique non-runaway P-symmetric solution at that given positive energy. Numerical 
analysis^ shows that at some critical energy this unique solution splits into three 
solutions, from which two are not P-symmetric themselves, but are P-conjugated 
to each other. 

For the 3D WF the only known result is the linear stability analysis of the 
Schonberg-Schild circular orbitaiS, exhibiting an infinite number of unstable solu- 
tions for the characteristic equation. This paradoxical linear instability of circular 
orbits is a warning that the first relativistic correction can not be the generic un- 
folding of this complex dynamics considered as a perturbation of the integrable 
Coulombian system. The question whether these unstable solutions are present in 



Wheeler- Feynman electrodynamics 3 



the exact theory can be investigated numerically, which is a use for our integrators. 

In the following we discuss three numerical methods for the solution of 3D 
WF. The first one, originally proposed in EH, involves solving the equations with 
respect to the most advanced velocities and accelerations, converting the equation 
of motion to a double-delay form. As a result, the evolution becomes completely 
and quite naturally determined by the past trajectory, allowing an easy numerical 
implementation. However, in section 2.1 we show that this numerical scheme is 
largely unstable and can not be used in practice. We discuss the reasons of this 
instability, which are closely related to the existence of unstable solutions of the 
characteristic equationlill In section 2.2 we describe a method, based on usage 
of special parametrization of world lines (light-ladder parametrization □), which 
implements the direct integration of the equations when the integration step is equal 
to one step of light ladder. The method is applicable at low velocities v/c ~ 1CP 2 
and is stable at least for integration times up to N ~ 10 revolutions of the particles. 
In section 2.3 we consider an iterative method, originally proposed in 111 to solve the 
1-dimensional 2-body WF problem, which after certain modifications can be applied 
to the 3-dimensional case. This method works with shorter intervals of integration 
than the light-ladder and it is able to resolve the structure of solutions for intervals 
smaller than one light ladder step and also converges at high energies. The results 
obtained using the methods are presented in section 3. 

The mathematical structure of the equations of motion of 3D WF is described in 
Appendix A. Poincare's invariance of the Fokker Lagrangian implies the associated 
non-local .Moether's constants of motion as integrals over a segment of the past 
historyOuta. These conserved integrals provide a useful check for the numerical work 
and we discuss them as well in Appendix A. Further appendices provide necessary 
details about numerical schemes. 

2. Computational methods 

2.1. Solving the equations with respect to the advanced variables 

The most intuitively attractive method, which is based on solving the equa- 
tions of motion for the most advanced velocities and accelerations, turns out to be 
numerically unstable. The reason for this instability can be easily understood con- 
sidering the simpler 1-dimensional version of the problem with repulsive potential. 
We henceforth use the letter x to denote the cartesian vector of position for particle 
1 and the letter y to denote the position of particle 2. The velocity and acceleration 
evaluated at the present time of particle 1 are written as v x and a x respectively. 
An upper plus (or minus) above quantities signifies that the quantity is evaluated 
at the future (or past) point of the light cone starting from the present position 
of other particle. The equations of motion for 1-dimensional repulsive motion have 
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the simple form 

a x = 1(1 - ^)3/2 fI±V . {y - _ x) -2 + i_< . ( + _ x) -2\ 

\1 — Vy 1 + Vy J 

in units e = m = c = 1. The equations for the other particle are obtained by 
exchanging x <-> y in the above formula. The idea is to solve the above equation of 
motion of particle x for the most advanced velocity Vy of particle y as 
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The equation for v+ is again obtained by interchanging x «-> y. It is seen that we 
have a functional differential equation with delayed argument only, as first pointed 
out in0. As far as initial conditions go, the general theory on delay equations 
tells us that we need to provide an initial C 2 function describing the position of 
particle x in the past, as well as the information on particle y needed is also to be 
provided over twice the retardation lag seen by particle y. 

The first difficulty with such setting of the problem is that the expression for 
v£ includes the derivative a x — dv x /dt x evaluated in the retarded time. In the 
numerical solution this derivative should be estimated from known values v x (t x ) 
on the integration grid using an appropriate finite difference scheme. The common 
feature of these schemes is that they contain a factor h^ 1 , where h is the integration 
step. Supposing that points v x (t x ) on the grid have a computational error of order 
e, we will have an error in the derivative dv x /dt x of order eh^ 1 . Due to the func- 
tional relationship u+ = f(dv x /dt x ), Vy will have the error eh^ 1 . We see that the 
numerical scheme amplifies the computational error by a factor of h^ 1 after each 
step of light ladder. The integration step for the given equation should be small, 
typically h < 10~ 2 , and as a result, the scheme works only for 2-3 light ladder steps 
and then explodes. 

The second difficulty concerns the possibility of continuation of the trajectory. It 
turns out that the advanced velocity given by the above expression is not guaranteed 
to be bounded by the velocity of light. Even in the case where numerical integration 
would be possible, for some initial data, at a certain point Vy can exceed the velocity 
of light and the solution cannot be continued to the next light ladder step (where 
this v y > 1 should enter under the square root (1 — v 2 ) 1 / 2 ). At least for symmetric 
1-dimensional orbits Driver's theoremEJ says that the set of infinitely continuable 
solutions is finite-dimensional, implying that for most initial functions the solution 
can not be continued for all times, and only a thin finite-dimensional subset of initial 
functions produces the good solutions. Ideally, in the given approach one should 
find special rare initial data, for which the solution can be continued for a large 
time, and this is easily seen to be a formidable numerical task. 

These same difficulties plague the 2- and 3-dimensional cases, with the extra 
complication that the equations of motion include the advanced acceleration in a 
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degenerate way. An algebraic constraint appears and we solve it in Appendix B. 
As a result of this constraint, one component of the advanced velocity can be found 
from the previous acceleration (and other variables) by a functional relation of the 
form v + = /(a,...). This property again makes the resulting numerical scheme 
unstable, like in the 1-dimensional case. 

2.2. Direct integration scheme (low energy) 

In the limit of small velocities a step of light ladder is small and can be taken 
as an integration step. The following computational scheme can work in this limit: 




Fig.l. Direct integration scheme. 
Suppose that coordinates and velocities are known in points Xk,k < n; yu,k < 
n + 1 (shown by dots on fig.[l]a), and accelerations are known (stored from past 
integration) in points x k , k < n — 1 ; yu,k < n (shown by circles) . We then go 
through the following steps: 

51. Because the step of light ladder is small, we are able to compute the accelera- 
tion in point y n +i as the left derivative, using a finite difference scheme, described 
further. This computed derivative is shown by brackets in fig.[l]b. 

52. As a result, we are able to find the acceleration in point x n , using the equations 
of motion. This is shown in fig.[j]b by a circle around point x n . 

53. We are then able to perform one integration step, applying Stunner's method 
(described further), and find coordinates and velocities in point x n+ \. 

Next we apply these three steps for x <-» y interchanged. As a result, we have 
fig.[l]c with the same pattern of data distribution, as in figjlja, only with shifted 
index n — ► n + 1. 

Note: during this process we actually have two definitions of acceleration in each 
point: computed by current and past velocities as left derivative (shown by brack- 
ets), and computed from the equations of motion (shown by circles). These two 
accelerations should coincide, and their relative difference |Aa|/|a| can be used to 
check the consistency of the method (see Appendix C). 



Sturmer's integration method for the equation x — f(t,x) is described in 
@,p.20. This nth order formula can be obtained from the following Taylor expan- 
sions: 



<7fc = hf(t k ,x k ), 
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g fc _i = q k - x k h 2 + \xfh 3 - ... + o(h n ), (0.1) 

q k -2 = qk - 2x k h 2 + 2xfh 3 - ... + o{h n ), 

q k -n+l = ?fc - (n - l)x k h 2 + &=pl x W h 3 - ... + o(h n ), 

by solving these (n — 1) linear equations for the (n — 1) unknowns x^h p and 
substituting them into the Taylor series for x k+ i 

x k+1 = x k + q k + \x k h 2 + ^x^h 3 + ... + o(h n ). 

Particularly, for n — 8 we have 

x k+1 = x k + q k + \&q k -i + j|AV._2 + |A 3 gfe-3 + 2 ^A 4 q k ^ i 

_i_ 95 A 5 i 19087 a 6„ , 5257 a7„ , „fu8\ In n\ 

+ 288 A q k -5 + 60480 A + mso A q k - 7 + o{h ), (0.2) 

where A p q k are defined recurrently as 

A 9fe _i = q k - q k -i, A p q k -i = A p ~ l q k - A p ~ 1 q k -i. 

The second derivative x k , found from (^l]) at n — 8, is given by the following 
backward differentiation formula: 

7 

x k = h- 2 J2l AP <lk- P + o(h 6 ) (0.3) 
p=i 

Notes: 

1. The advantage of Stiirmer's method is that the right hand sides of the equations 
do not need to be estimated in the intermediate points of the grid x k , as for example 
with the type of Runge-Kutta methods, thus the necessity of interpolation of data 
between integration points is avoided. 

2. Usage of lower order scheme reduces the precision of integration, but improves 
stability of the method. Particularly, the reduction of order n = 8 — * n = 4, 
correspondingly reduces the precision of conservation of Noether's integrals (see 
appendices), and increases the upper velocity, for which the method is convergent 
v = 0.006 -> v = 0.02. 

2.3. Iterative scheme (high energy) 

The following scheme, originally applied in El to solve the ID WF with repulsive 
interaction, is here improved and implemented to integrate the 3D WF in the highly 
relativistic limit. The method works as follows: 
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An initial trajectory is guessed for the two particles, 
which is called iteration number zero. We then use 
the equations of motion to evaluate the present force 
caused by the past data points x~ ,y~ , located on 
current trajectory, and future data points x + ,y + , 
taken on the trajectory of previous iteration. This 
acceleration is used as a force field to integrate the 
trajectory, producing the next iterated trajectory. 
Repeating such integration of trajectories for fixed 
data on the initial segment, we arrive (if the itera- 
tions converge) to a solution of the_WF problem for 
the initial data. As mentioned in 113, this scheme is 
not guaranteed to be convergent, and we developed 
special methods of stabilization (see Appendix C), 
that is able to keep convergence up to very high val- 
ues of the velocity. 



y x / 




Fig. 2. Iterative scheme. 



We start our iterations from the analytically known circular solutions^, fixing the 
initial segments of trajectories r £ [0, 0.5] to arcs of circular orbits through all 
iterations. Then, to obtain non-circular orbits, we apply short pulses of external 
force to both particles immediately after the segment of initial setting. The method 
allows to subdivide each light ladder step into a large number of integration steps 
of sufficiently small size to have the correct integration at high energies. Stunner's 
method, described in the previous section, can be used for integration. Special 
boundary conditions should be used at the end of integration interval, where no 
information is available about the advanced Lorentz force. These and further details 
about the method can be found in Appendix C. 

This iterative method requires storage of the complete integrated evolution from 
the current and previous iterations. Due to memory restrictions, only a short part 
of the trajectories (up to 10 revolutions) can be determined by this method. Nev- 
ertheless, the method can provide useful information about structure of solutions 
at high energy when only short term evolution is needed. 



3. Results 

The computational methods, described in the previous sections, are applicable 
to arbitrary masses and charges. However, the case of equal masses possesses ad- 
ditional symmetries, and here we present the results obtained for equal masses and 
opposite charges (positronium-like system). Following Andersen and von Baeyerllj 
we fix the unit system to c = 1, e x = — e v = 1, m x = m y = 2, so that reduced mass 
m = (m~ + rriy 1 )^ 1 = 1. We present the evolution of the system in the center- 
ofpmass frame (CMF), using Schild's definition of CMF origin, given by Eq.(2.13) 
irM 
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haw energy: Darwin's precession. At low energy the solution looks like very 
slowly precessing ellipse, shown on figj^a. Between each two consequential ellipses 
of this figure 500 intermediate ellipses are not shown. The figure presents two 
coincident curves: trajectory of particle x and trajectory of particle y, reflected 
about the origin in CMF, demonstrating the P-symmetry of the solution. 
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Fig. 3. Low-energy solutions: a) trajectories, b-d) related parameters. 



Related parameters (E, L, loq, uj p ): binding energy E = 2m — yP^, angular 
momentum in the CMF, angular frequencies of rotation and precession. For 10,000 
solutions generated by the methods described in Appendix C, we measure these 
four parameters, and analyze the 4-dimensional scatter plot. Figures |3p,d represent 
projections of the scatter plot onto planes (E,L) and (E,Wp/wo). The boundary 
curves on these figures correspond to the circular orbit limit, and are given by 
formulae Eqs.(3.4),(3.5) in E3 and Fig. 4, Eq.(5.1) in LL3, which for the case of equal 
masses and low velocities reduces to E = 1/(2L 2 ) = uo p /uo . 

Figure ||b represents a diagram (L,oj p /ojo), where the scatter plot is projected 
almost to a curve. It shows that the ratio ujp/ujQ depends mainly on L as uj p /ujq = 
1/(2L 2 ). This dependence was predicted by Darwin E3, who first studied the lowest 
order relativistic corrections to the Coulomb problem and found the described effect 
of precession. Using non-perturbative methods, we detect deviations from Darwin's 
law by an order 10 -10 , still resolvable by our measurement (which gives absolute pre- 
cision of u>p/u!o about 10~ 14 ). Inner image shows the difference Au> p /uio = 1/(2L 2 ) — 
ujp/ujQ versus L, for 4 narrow bands 2EL 2 = [0.98,1], [0.78,0.8], [0.58,0.6], [0.38,0.4], 
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following on the image in the order down-up. Thus, u) p /loo depends mainly on L 
and very weakly on E. 

The structure of the graphs shows that the scatter plot in 4-dimensional space 
(E, L,ujq,uj p ) is located on a 2-dimensional surface, like in the non-relativistic 
Coulombian problem, where (E, L) are the only parameters, defining the shape 
of solutions in the CMF. Another similarity with Coulombian problem is that the 
trajectories in CMF with high precision* belong to a plane, perpendicular to the 
vector of angular momentum L, even if initial data in laboratory frame are non- 
planar. In relativistic 3D WF problem we cannot prove planarity of orbits so easily 
as in Coulombian problem, and currently we have it only as a result of numerical 
experiment. 

High energy: Bifurcations. Figure || shows pieces of the high-energetic trajec- 
tories, found by method sec. 2. 3. Here we again superpose trajectories x and (— y) 
to test the P-symmetry of the solution. On the left image a small asymmetry is 
visible at the beginning of the trajectories, caused by the asymmetrical statement 
of the problem at the beginning (see Appendix C), which exponentially disappears 
in the inner regions of the trajectories. When the energy increases, the asymmetry 
penetrates deeper to inner parts of the trajectory, and at high energy the solution 
completely looses the symmetry. 

We define the measure of asymmetry as r\ — max|x(i) + y(t)\ in CMF, where 
the maximum is taken over a step of light ladder. Then we draw graphs rj 2 (E) 
for consequential light ladder steps - they are given by curves with points on the 
left graph of figj|. We see that at low energies the asymmetry tends to zero with 
increasing number of light ladder steps, while at high energies the asymmetry tends 
to a non-zero value. The limiting ij 2 is well described by a linear function of E 
(corresponding to a supercritical pitchfork-type bifurcation r/ ~ ■sfE — Eq), it is 
shown by upper straight line on the graph. Lower straight lines correspond to 
analogously constructed dependencies for smaller asymmetry of initial conditions. 
Right image presents the result of assembling all pictures to 3D graph. On this 
graph AL(E) — L(E) — Lq(E), Lo(E) is angular momentum for circular orbits. 
Line AB subdivides the diagram onto two parts: symmetric phase - on the left 
from AB the trajectories possess P-symmetry and the shape of solution is uniquely 
defined by two parameters (E, L) ; asymmetric phase - on the right from AB the 
symmetry is violated and given (E, L) correspond to a continuous set of solutions, 
differing by measure of asymmetry. 

Such a structure of the phase space was actually predicted by Andersen and 
von Baeyerta, from the linear stability analysis of the circular orbits in 2D WF. 

*Motion along L-axis consists of small oscillations with amplitude Az = 10~ 8 of classical radius 
or relative value Az/(XY-size of orbit) = 10 — 12 ; and slow Brownian motion with the velocity 
about v = 10 — 17 ; both motions have the same order as non-conservation of Noether's integrals, 
appearing due to computational errors. 
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This work shows that at velocities of about v — 0.95, new real eigenvalues of the 
characteristic equation appear, representing new degrees of freedom in the system. 
Corresponding linear mode, denoted by Andersen and von Baeyer as D^', violates 
the P-symmetry. In our non-perturbative approach we detect this effect in the 
region of energy E E [2.6,2.8], which for circular orbits corresponds to the region 
of velocities v € [0.937, 0.954], in good agreement with the predictions Actually, 
in point B, where the orbits are close to circular, we have E = 2.75, v = 0.950, i.e. 
exact coincidence. 
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Fig. 4. High-energy solutions: trajectories. 
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Fig. 5. High-energy solutions: related parameters. 



Conclusion. Two numerical methods are developed to solve 3D WF problem in 
the shallow and high energy ranges. The methods go further than pre-existing 
onestEftj, based on first order perturbation theory. Our non-perturbative analysis 
reproduces part of the pre-existing results with the exception of the divergent modes 
predicted in 0, whose exponential increase violates the limitations of the linear 
perturbation theoryEl. The question about their existence in the complete theory is 
open. Our numerical methods do not reproduce such solutions as far as the unforced 
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equations of motion are considered. On the other hand, some traces of divergent 
modes are present in the solution of forced equations, describing the reaction to an 
external perturbation. The results of this analysis can be found in E3 and will be 
published elsewhere. 
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Appendix A: equations of motion and conserved quantities 

1. The equations of motion for the 3D WF have the following form: 



a 



(±) - ^v^i { m [ i - + (t - * ) ) . (0.4) 



± I I \^x "X, I ± 
r x0 J \ r x0 

r- — I <r 1 1 1 W 



ft=y ± -x, r± =±\rt\, a x = ±(ai +) + 4 _) ) + vl{F ex t - {F ex tVx)v x ). 

mix 

Here a x denotes the acceleration of particle x; variables E^r are advanced and retarded components 
of electric field, created by particle y at the location of particle x; variables y^- , Vy~ , a y stand 
correspondingly for position, velocity and acceleration of particle y at time moments, where the 
world line of particle y is intersected by the light cone with origin in x. Equations of motion for 
particle y are given by replacement (x <-> y) in these expressions. The term F ex t corresponds to 
the force of external perturbation, used to produce perturbed initial conditions starting from the 
exact circular solution. For the external force we choose a cap-like form F ex t{t) = -F a exp[(l — 
(Ai/(t — to)) 2 ) -1 ], if \t — to I < Ai; otherwise 0. This function is infinitely differentiable and has 
a finite support. 

2. In our approach the equations of motion are integrated in the light ladder parametrization, 
which was introduced in [j and modified by us in the following way: 

On the initial segments of the world lines, shown in bold on fig.g, we . v / >yi-^— 

select arbitrary monotonous parametrization r S [0,0.5], e.g. taking r ~ 
as a linear function of the time components xo,yo- We then consider 
light rays, which consequentially reflect between the world lines, and 
distribute the parametrization along these rays, adding r — > r + 0.5 
at the points of reflection! In this parametrization x(t) and y(r±0.5) 
are related by light rays, so that deviation of the argument in the 
equations of motion becomes constant. As a result, while determining 
the intersection of the light cone with a world line, we are taking the 
data directly from one of the past integration points, not solving the 
equation of intersection and not performing high order interpolations Fig. 6. Light ladder 
of data between integration points. parametrization. 

Using the light ladder parametrization, we consider the evolution of the variables x(t),j/(t), 
xo(r),yo(r) reformulating the equations of motion in terms of x = xov x , v x = xoa x , and the anal- 
ogous pair for the other particle y. We substitute for xo a formula xo = y^ ( |r^T | + {r^Vy )) / {\r^\ + 
(r^v x )), which follows from differentiation of the light cone condition r~^dr~^ = and expresses 
xo in terms of the past variables (y~ , Vy , i/^ ) and current variables (x , v x ) . For a x we substitute 




^This procedure gives a C°° -parametrization under special conditions in the vicinity of the point 
r = 0.5, which are satisfied particularly when the initial segments of the trajectories are set to 
circular orbits. 



12 I.N.Nikitin and J.De Luca 



the equations of motion given above (variable a,y , participating in r.h.s. of the equations, is taken 
from the stored past history, while Sy is taken in the iterative scheme of sec. 2. 3 from the trajectory 
of the previous iteration and in the direct integration scheme of sec. 2. 2 as a y = v y /yo, where v y 
is defined by backward differentiation formula). After these substitutions, we have a closed sys- 
tem of differential-difference equations for x(t), y(r). Variables xq(t), yo{r) do not participate in 
the equations, and can be determined by known solution of differential-difference equations using 
recurrent formulae xq = y^ + \y~ — x\, yo = Xq + \x~ — y\. 

3. The conserved quantities — Noether's integrals of motion, given by Eq. (2.9-10) in 0, after 
partial integration can be written in the following form: 



+~ A+(y-) - \{A-{x)A+{y-))(x - y~) 
-G^\x) + G^\y-), P^ = (x^y), 




= V / dr^F+^x), G<»>(0 = (x «-» y). 

2 J Q Fig. 7. Definition of 

Noether's integrals. 

Here j4 (j (^) everywhere means vector-potential felt in point (£). For example, A^(x) are the 
advanced and the retarded components of Lienard-Wierchert potential, created by particle y } felt 
in point x. The explicit formula for the potentials: (x) = ±e y yjj; /(r^y^). Analogously, F^ v {x) 
is the field tensor, created by particle y, felt in point x, given by formula F^ l/ (x) = 8A^(x)/dx^ , 
where square brackets denote antisymmetrization of indices. The above defined G-functions are 
equal to the integrated advanced Lorentz forces, acting on particles x and y. Using the fact that 
a sum of advanced and retarded Lorentz forces after integration give a local expression (change of 
momentum), G-terms can be rewritten to various equivalent forms, however it is not posible to 
reduce them completely to a local form (integrals cannot be eliminated). In our numerical scheme 
the G-integrals are computed in parallel with the main integration processes, by the same high 
order scheme. The given expression for Noether's integrals is graphically displayed by a diagram 
fig.pl, where the circles mean m^y-terms, solid lines are the segments of integration in G-terms, 
and arrows the denote potentials (directed from source to destination). 

This representation has an important difference from the Noether's integrals of the non- 
relativistic case (Galilean): The total momentum is composed of two contributions PS and 

Pa , where P^ depends on the position of the triple of points (y~xy + ) and P^ depends on 
the triple (x~yx + ). These triples are independent, as a result, both contributions are conserved 
separately: Pjf^ = Const, P^ = Const. It's easy to check that differentiation of p^^' 1 gives 
separate equations of motion for particles x and y. 

The expression for the angular momentum tensor also consists of two separately conserved 
contributions L M „ = ife + L^J : 

= + yff^On - G $(*) + g&V). J$ = (* « v), 

V„ W = m x ^= + ^-A-(x) + \{A-(x)A+(y-))yZ, vt) v) = (x « y), 



G<$ (?) = y f dr x x p x {li F+ v] (x), G$ (?) = (*«»). 
Jo 

The above defined Noether's integrals are used to define the centet|-of-mass frame (CMF). 
Particularly, defines the world line of CMF origin by Eq.(2.13) in £3: c M (A) = L pv P v /P 2 + 
XPfj,. In the coordinate system where the time axis is directed along this line and the origin is 
located on this line with c; = 0, three components of the angular momentum tensor form the 
angular momentum vector: L = (L23, L31, L12) and the other three components vanish: Lq; = 0. 
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Check of this property, as well as conservation of Noether's integrals were used to control the 
precision of numerical analysis (see Appendix C). 

Appendix B: Solving the equations of motion with respect to the most 
advanced variables 



1. From equation ( |),4| ) we see that the advanced contribution to the acceleration a x is linear on 

(1 - -i^) 1/2 ((l - nv x )8ij + (m - v x i)v x j), 



the electric field E. 
„(+) 



AijE+j, Aij(y+;x,v x ) 



e x 



where n 



"I, 



An analogous expression can be written for a 



(-) 



It is 



possible to prove that matrix A is non-degenerate, and that E x can be unambiguously found as 



Taking a 



(+) 



2d x 



we finally express E x as a function of the "future" 



variable y+\ present variables x,v x ,a x and past variables y 



1. According to (0.4), the electric field E x depends on the position of the source y + , its velocity 
Vy and its acceleration a,y , and has the following cartesian components: 



ZfT — ( n ~ v t)i + ((! - n v i)5ii ~ (n - a+. 



Matrix B is degenerate, i.e. riiBij = Bij(n — Vy)j = 0. In order to solve this degenerate linear 
system with respect to dy , we have to comply with the following: 

a) condition of consistency (constraint found multiplying the system by rn): 



+ N2 



(1 - nvy~) 



+ \2 



(0.5) 



b) When this condition is satisfied, the solution can be written as 



l(l-nv+fE+ + \(n-v+) 



(0.6) 



where A is a so far arbitrary parameter. In the 2-dimensional case, equation (1 — {Vy) 2 )/{1 — 
nVy~) 2 = C with C > defines an ellipses lying inside the circle (vy)" 2 < 1; and with C < - 
ellipses (—1 < C < 0), parabola (C = —1) and hyperbolas (C < —1), lying outside this circle. 
The physical region corresponds to C > 0, which means < 1. The consequence of the above 
is that it is impossible to solve uniquely for the advanced accelerations a,y from the equations 
of m otion Instead, we find a one-parameterr-4egenerate solution plus one constraint: equation 
(0.5) imposed on the velocity. In the followingCJ we choose a convenient set of variables to resolve 
the constraint and make the system availabkrfor explicit integration, then perform the integration 
using an adaptive delay equations integratorE3. In the result we find instability described in section 
2.1. 

2D case 
i'y plane 



Fig. 8. Algebraic constraint. 
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Appendix C: Details on the statement of the numerical experiments 

In the direct integration scheme the initial segments of trajectories were set to circular orbits 
This setting made in one half of one ladder step r £ [0,0.5] is sufficient to define the further 
motion. However, using Stiirmer's method of integration, we need also to initialize the finite 
differences APq^. Due to this reason we are setting initial segments of trajectories to circular 
orbits on longer intervals (up to 8 integration steps) . After this setting the integrator recovers the 
correct circular motion. To consider non-circular orbits, at a short time after the initial setting we 
apply to one particle a local external force F ex t- As a result of this action, the system acquires 
non-zero total momentum and starts to move away from the origin. We then transform back to 
CMF and observe the evolution there. 



Fig. 9. Initial segment of trajectory: 
a) in direct integration scheme; b) in 
iterative scheme. 




Usin g St iirmer's formulae of 8th order (0.2) and the corresponding backward differentiation 
formula (D.3), we have the acceleration, estimated with a precision of o(/i 6 ). The right hand side of 
the differential equations f(f^,X},, ...) are dependent on the computed acceleration, and it also has 
precision o(h 6 ). As a result, in the 8th order Stiirmer's formula: x^+i = %k + hfk + ... + o(h s ) the 
last term (A 7 gj,_7 ~ h s ) becomes ineffective, actually the method has 7th order. However, the 
term A 7 qj._7 should be in any case computed to find the acceleration and can be preserved in this 
formula. Note also, that the acceleration terms in the Lorentz force are suppressed by additional 
small factor e 2 /mr = ro/r (ratio of classical radius to particles separation) so that true precision 
of the method is between o(h s ) and o((ro/r)h 7 ). 

Using the grid of the light ladder parameter r = n £ Z, we still have convergent scheme, until 
the higher order differences A k q n _k are small. Actually, the role of parameter h here is played 
by retardation angle 8 ~ v/c. Near circular orbits we have an estimation ro/r ~ (v/c) 2 , so that 
for v/c ~ 10~ 2 the application of 8th order scheme gives the error of order 10~ 16 , on the level of 
machine double precision. 

The integration process can be considerably accelerated by means of caching of previously 
computed data. For this purpose we store a short buffer of past data (8 integration steps), 
including all variables, whose evolution is tracked, as well as the computed right hand sides of the 
equations of motion q^, their finite differences A p <jfc and some auxiliary variables (such as r x y , 
which due to relations /^(t) = —Ty{r ± 0.5) enter in computation several times). The content of 
the buffer is periodically written to a file to allow restarting of the program from previous states. 
To transform the system to CMF we perform a Lorentz transformation of all variables in the buffer 
(some of these variables, e.g. accelerations, have non-linear transformation laws). 

Control of the precision was performed by comparing two definitions of acceleration (variable 
AS in sec. 2. 2) and testing the conservation of Noether's integrals (Appendix A). The relative 
difference [Aa[/[a| was < 3 ■ 10 -10 in the region of free motion (in the region, where the external 
perturbation is applied, this value was about 0.07, i.e. the method has worse precision in the region 
of perturbation, but further free motion is integrated precisely). As already mentioned, Noether's 

(x) (v) 

integrals of Appendix A separates in two conserved quantities P M = +Pjf . Their conservation 
was satisfied in our numerical scheme (order n = 8) with a precision of |AJq |/Po < 10~ xl , 
\&p(*-y)\/P Q < 10~ 13 . For total values in CMF we had \P\/Pq < 10" 14 . For the angular 
momentum tensor we had found that its components have a cumulative numerical error, and for 
these values in the CMF we have an estimation \L i\/\L$\, \A£\/\L\ < 10~ 14 *num.of revolutions, 
valid up to the maximal number of revolutions = 10 6 we considered. 

The direct integration scheme is applicable upto v = 0.006 (n = 8). For greater velocities 
the solutions exhibit oscillations with time size comparable with the integration step, and rapidly 
increasing amplitude. The appearance of these oscillations is followed by violation of the conser- 
vation laws. This effect is caused by the increase of the light ladder step to a critical value, when 
it cannot be used anymore as a step of integration. The upper limit of v creates on figbb the right 
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limit of E (high velocity almost circular orbits) and lower limit of L (in this limit the trajectories 
are thin ellipses, and the particles come close to each other, reaching the upper limit of v). 

In the considered range of velocities the precession of the orbits is so slow, that the measure- 
ment of its angular frequency u> p =(angle of precession)/ (time) requires special efforts. For this 
purpose we find intersections of the orbits with the circles r = Const and measure the difference of 
the angles Aa; for consecutive intersection points. For better precision we interpolate the solution 
between the integration points using the same high order polynomial of the integration method 
and find the intersection using dichotomy. Then we use large statistics of Acti to find the average 
Aa, which is related to the frequences as Aa = 2tt(ui p /ujq). Its mean-square error consists typi- 
cally of about 10 — 13 radians, leading to a relative error of (a)p/wo) of about 10 — 9 , sufficient for 
our purposes. 



Fig. 10. Measurement of precession rate. 




Construction of the scatter plot fig.g: The four main input parameters of the program are 
the radius of the initial circular orbit and the 3 components of the external perturbation force. 
These were varied in a certain 4-dimensional region and 10,000 points w ere generated in the run. 
In 25% of cases the system was ionized (all these cases correspond to V P 2 > 2m, i.e. when the 
external force gives the system enough energy for ionization). In 0.8% of the cases the particles 
were too close to each other and moved too fast and as a result, the algorithm lost applicability 
and diverged. These cases were rejected at earlier stages of the integration and did not consume 
much computation time. The remaining cases are displayed on fig.pl 

To prove that the scatter plot occupies a 2-dimensional surface in the 4-dimensional space, 
we consider its sections by hypersurfaces EL 2 = Const. The projection of these sections on the 
plane (L, Auj p /u>o) is shown on the inner image of fig j^b; projections to the planes (E, Alu p /u>o) 
and (luq, Ali) p /u)q) look similarly. As a result, we can conclude that the sections are 1-dimensional 
curves in a 4-dimensional space, which span a 2-dimcnsional surface with a change of the section 
parameter. Two other approaches to measure the dimension of the scatter plot are described in 
r% (1) computation of the Jacoby matrix for the mapping from the space of input parameters 
to the space of measured characteristics, and estimating its rank; (2) direct visualization of this 
4-dimensional object by projection to the 3-dimensional space with color encoding of the 4th 
coordinate. 

2. In the iterative scheme the initial segments of the trajectories r £ [0,0.5] were set to circular 
orbits. To obtain non-circular orbits, we apply to both particles the local external forces —F ex t and 
Fext + AFext- Such setting is convenient to control separately the orbits deviation from circular 
and asymmetry of initial conditions (the first one is controlled by F ex t, the second by AFext). 

Note: There is still another type of the initial condition, where the data on initial segment are 
directly set to non-circular orbits. In principle, these two types of initial conditions are equivalent: 
the equations of motion are generally not satisfied on the initial segment, where the trajectories are 
set to arbitrary shapes, and non-zero right hand side of the equations in this segment play the role 
of an effective external force. However, concrete implementation of the non-circular initial condi- 
tion creates certain problems in our scheme. The data on the initial segment should necessarily be 
presented in light ladder parametrization, for which the analytical representation is known only for 
the case of circular orbits. Moreover, the correct parametrization on the whole trajectory is recov- 
ered by further integration process from the initial segment, and its C n -smoothness is guaranteed 
only if special conditions are satisfied (left and right derivatives upto nth order should coincide 
in a point r = 0.5). For non-circular initial orbits the satisfaction of this property is problem- 
atic, while for circular orbits a global C°°-smooth parametrization can be easily constructed, and 
infinitcly-diffcrentiablc external forces, used to control the shape of solution, preserve smoothness. 

For convergence of the method the total number of light ladder steps Nls and the number of 
integration steps per light ladder step Nps should be set to Nps > 10 2 , Nls < 10 at intermediate 
energies E ~ 1 and to Nps > 10 3 , Nls < 5 at high energies E ~ 3. Our stop-criterion of iteration 
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was J2 i ((%ri(n)-x n -i(T i )) 2 + (y n (Ti)-y n -i(Ti)) 2 )/(Nps*Nls) < e, where (x n ,y n ) represent the 
current iteration, (af n _i,y«_i) - previous iteration, the sum is taken over all integration points, 
and e = 10" 12 . 

In the described simplest form the method can be applied up to velocities v = 0.8 1), then 

a stabilization is needed to keep it convergent. A stabilization procedure, proposed in E3, suggests 
to combine the acceleration computed by the equations of motion with the acceleration stored from 
the previous iteration to a n (t) = 0.1a eq (t) + 0.9a n — l(i), and integrate it to the trajectories of the 
current iteration. We have found that this type of stabilization, being applied to our problem 
at high energies, makes the method indifferently stable: it has no divergence but cannot reach 
the solution either. We have used another approach, based on the following predictor-corrector 
algorithm: We track the solution as a function of one control parameter /i (which, for example, 
can be equal to initial velocity of the particles). The solutions found for the n previous values of 
fj, should be stored (n old copies of coordinate, velocity and acceleration arrays should be kept; 
while other variables, entering in the caching buffer, should not be copied). The dependence of the 
solution on /i at each integration point is approximated by a polynomial of fcth order using a least 
square method (in our implementation n = 10, k = 2). This polynomial is used as a starting point 
for the iterations at a new value of the control parameter /i + Afi. If the solution has a smooth 
dependence on the control parameter, this method produces a starting point placed much closer 
to the result then other types of initial guess. This fact considerably improves the convergence. 
The step Afi is chosen adaptively: decreased twice if convergence is lost and increased twice if the 
solution has beenrfound, automatically keeping this value in the optimal region. This approach, 
earlier applied in EJ for the solution of ID WF, accelerates the recovery of the solution along the 
^t-axis by a factor of thousand, and keeps the method stable up to the velocity v = 0.98. 

Other parameters, influencing the shape of solution, such as the components of the external 
force, in this method should be set to the given smooth functions of fi, which should change 
from zero to (for example) certain constant values. Starting the procedure for circular orbits at 
the velocity v = 0.8, where the simplest scheme is convergent, we track one-dimensional families 
of solutions, corresponding to the given functions, and finally cover the whole space of input 
parameters in the region of convergence of the method. On fig.0 the right limit in E, bottom 
limit in AL and upper limit in r] 2 are the limits of convergence of the method respectively at 
high energy, strong perturbation force F ex t and large asymmetry AF eX f The upper limit in AL, 
close to the circular orbit, is related with a problem of another kind: Such solutions have very 
good convergence, but do not exhibit linear dependence in graphs rj 2 (E), so that we cannot clearly 
separate symmetric and asymmetric phases by the described method. 

An extra feature of the iterative scheme is the existence of artificial boundary effects, appearing 
near the end of the data arrays. Because at the end of the integration interval the future evolution is 
not known, the advanced Lorentz force should be omitted in this region. For the retarded Lorentz 
force we have used two possibilities: (a) single retarded Lorentz force or (b) double retarded 
Lorentz force is taken on the last step of the light ladder. Absence of the advanced force and 
various settings for retarded force lead to boundary effects, propagating to the inner regions of the 
integration interval. The amplitude of these effects is exponentially decreasing with the depth of 
penetration, so that the solution at several light ladder steps before the end of the integration is 
no more sensitive to boundary effects. It has also been found that boundary conditions of type (a) 
lead to stronger boundary effects, but support better convergence at higher energies than condition 
(b). Figjll| shows the solutions, corresponding to boundary conditions of both types for value of 
velocity v = 0.86. 



Fig.ll. Two types of boundary ef- 
fects in iterative scheme. 

a) b) 

The described modification of the Lorentz force leads to non-conservation of the integrals of motion 
in the last step of light ladder. At the beginning of the trajectory in those regions, where the 
external forces are applied (and where the initial data are set — as we show above, this setting 
leads to the appearance of an effective external force) Noether's integrals are also not conserved. 
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On the remaining parts of the trajectory the conservation laws are valid, namely these parts are 
shown in the images fig.UJ. The conservation laws in these parts were satisfied with precision 
\P\/Pq < lCT 10 , |AL|/|L[ < 10~ 8 , |L i|/|f| < at intermediate energies E ~ 1, slowly 

increasing upto |P|/Po, |AL|/|L| < 10" 6 , |Loi|/|£| < 10 ~ 4 at hi ~ h 

energies E ~ 3. 

3. Computational time: most extensive computations are required for the construction of the 
scatter plot figfl (24 hours) and 3D graph fig.0 (120 hours). The computation was performed 
during one night, parallelly on 12 processors (300MHz MIPS R12000) of the GMD computer 
SGI/Onyx2. 
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